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Abstract — Design of Random Modulation Pre-Integration sys- 
tems based on the restricted-isometry property may be subopti- 
mal when the energy of the signals to be acquired is not evenly 
distributed, i.e. when they are both sparse and localized. 

To counter this, we introduce an additional design criterion, 
that we call rakeness, accounting for the amount of energy that 
the measurements capture from the signal to be acquired. 

Hence, for localized signals a proper system tuning increases 
the rakeness as well as the average SNR of the samples used in 
its reconstruction. Yet, maximizing average SNR may go against 
the need of capturing all the components that are potentially 
non-zero in a sparse signal, i.e., against the restricted isometry 
requirement ensuring reconstructability. 

What we propose is to administer the trade-off between rake- 
ness and restricted isometry in a statistical way by laying down 
an optimization problem. The solution of such an optimization 
problem is the statistic of the process generating the random 
waveforms onto which the signal is projected to obtain the 
measurements. 

The formal definition of such a problems is given as well as 
its solution for signals that are either localized in frequency or 
in more generic domain. 

Sample applications, to ECG signals and small images of 
printed letters and numbers, show that rakeness-based design 
leads to non-negligible improvements in both cases. 



I. Introduction 

This paper is about the application of some recently devel- 
oped signal-processing techniques to the sensing of physical 
quantities, i.e., to their conversion into a sequence of samples 
that can be processed by an electronic system for the most 
diverse purposes. 

Conventional approaches to this are based on the celebrated 
Shannon-Nyquist theorem stating that the sampling rate must 
be at least twice the highest frequency in the band of the signal 
(the so-called Nyquist frequency). This principle is the basis of 
almost all methods of acquisition used in nowadays audio and 
video consumer devices, in the processing of medical images, 
in the operation of radio receivers, etc; 

Compressed Sensing (CS) is a recently introduced paradigm 
for the acquisition/sampling of signals that violates the 
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Shannon-Nyquist theorem providing that additional (actually 
surprisingly broad) assumptions can be made. 

A bird's eye view of CS shows that it is based on two 
general concepts: sparsity, which materializes the needed 
additional assumption, and incoherence between coordinate 
systems. 

Sparsity expresses the idea that the information content of 
a signal can be much less than what is suggested by his 
bandwidth, or, for a discrete-time signal, that the number of 
its true degrees of freedom may be much smaller than its time 
length. Actually, many natural signals are sparse in the sense 
that they have a very compact representation when expressed 
with respect to a suitable reference system and are therefore 
susceptible to CS. 

Incoherence extends the concept of duality between time 
and frequency. It is used to formalize the fact that when two 
domains are incoherent, objects that have a small representa- 
tion in the first of them spread their energy over a wide support 
when seen from the point of view of the other domain. 

It is evident that the first domain is best one when it comes 
to express and characterize the signal while the second is to 
be preferred for sensing operations since even few scattered 
measurements have a chance of capturing the signal energy. 
This is exactly what happens, for example, if we want to 
acquire a sinusoidal profile of unknown frequency. Since such 
a signal is extremely sparse in the frequency domain, the only 
two non-zero components of its spectral profile are incredibly 
effective in representing it. Yet, nobody would probe the 
frequency axis at few frequencies with the hope of coming 
across the one at which the signal is present and thus being 
able to recover the amplitude and phase. Instead, we know very 
well that only few samples in the time domain are enough to 
capture all the signal features. 

Generalizing all this, CS architecture analyzes the target 
sparse signal by taking few measurements in the domain in 
which the energy is widespread and thus easy to collect. If this 
is done properly, the resulting samples can be subsequently 
processed by algorithmic means to reconstruct the small 
representation in the domain in which the signal is sparse. 

The theoretical and practical machinery needed to do this 
in realistic conditions is being rapidly developed to arrive 
at acquisition mechanisms that can be labeled as Analog-to- 
information (AI) converters \\\. In fact, once that the proper 
domain has been found in the form of a waveform basis along 
which the signals can be expressed as a linear combination 
with a small number of non-zero coefficients, the actual 
information being carried by the signals will be found in the 
positions and the magnitudes of those coefficients |2][3| that 
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are the true target of the CS procedures. 

II. CS FOR LOCALIZED SIGNALS 

The leveraging on sparsity has been recently paired |4||5] 
with another technique widely used by engineers to spot 
information content in signals, i.e. the uneven distribution of 
average energy along properly defined bases (that, in general, 
are different from those for which sparsity can be identified) 


In the following, we will indicate such an uneven distribu- 
tion of energy with the term localization and we will observe 
that, in general, it provides a different a-priori information with 
respect to sparsity. As it is consequently naturally to expect, we 
will be able to show that these signal features allows improved 
sensing operations. 

The key assumptions under which this may happen are that 
(i) measurements are taken by projecting the signal onto a 
proper set of waveforms whose cardinality is smaller than the 
dimensionality of the signal, and (ii) the overall effect of dis- 
turbances in the sensing process (thermal noise, quantization 
errors, etc.) can be modeled as a projection-independent error. 

When (i) and (ii) hold, a noise-tolerant reconstruction of the 
sparse signal from a number of measurements that is smaller 
than its dimensionality is commonly achieved by designing 
the projection operator so that it is a restricted isometry (RI), 
i.e. it approximately preserves the length of the sparse signal 
to which it is applied so that the ratio between the norm of 
such signal and that of its projection falls within an interval 
[Vl - S, VT+~S] where the RI constant < 5 < 1 should be 
as small as possible J2) Q. 

Roughly speaking this means that, if the measurements 
come from a RI, the original signal energy is not lost in the 
projection and, when acquisition error is added, the signal- 
to-noise ratio (SNR) of the samples remains high enough to 
perform reconstruction. 

This approach and its pairing with localization can be 
intuitively explained with reference to a simplified, low- 
dimensionality setting in which the signal to acquire a has 
three components (a ,ai,a 2 ) and is sparse since only one 
of its components is non-vanishing in each realization. More 
formally we may assume that aj ^ with probability pj 
for j = 0,1,2 and, obviously, for po + pi + p-2 = 1. 
Furthermore, when it is non-zero, the j-th component of the 
signal is a realization of a random variable with variance cr| 
for j =0,1, 2. 

It is worth stressing that the possibility of po^o ^ Pi (J i 7^ 
P2&2 implies that localization and sparsity are two separate 
concepts. In fact, though a is sparse by construction, its 
average energy concentrates on the axis whose associated Pj<Tj 
is larger and this concentration depends on the unbalance 
between the probability-variances products. 

Since a is sparse, it can be reconstructed by measuring its 
projection on a two-dimensional plane. To define it, refer to 
Figure [T]-(a) and note that the generic projection plane passing 
through the origin defines an angle 9j € [0, 7r[ with each axis 

'A typical example is the class of band-pass signals, which are localized 
in the frequency domain, i.e., with respect to the Fourier basis. 




Fig. 1. A simple CS task using a projection plane designed by considering 
only the restricted isometry property (a). A graphical evaluation of the 
corresponding restricted isometry constant (b). 
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Fig. 2. A simple CS task using an optimized projection plane designed by 
merging rakeness and restricted isometry (a). A graphical evaluation of the 
corresponding restricted isometry constant (b). 



ctj for j = 0, 1, 2. These angles are such that Y^j=o sm2 ($j) = 
1. 

Any set of angles 8j ^ w/2 for j = 0, 1,2 is a feasible 
choice. This is shown in Figure [T[(b) that reports a unit radius 
circle on the projecting plane, along with the projections of 
unit-length segments centered in the origin and aligned with 
each of the three axis. 

As long as the three projections have no point in common 
but the origin (so that, in general, the projections of the 
coordinate axes are distinct straight lines) the retrieval of the 
original signal in noiseless conditions can be ensured without 
complicated algorithms. 

When noise comes into play, classical CS theory looks for 
planes corresponding to projection operators that are good 
RI. To do so, note that the ratio between the length of a 
segment aligned with the axis a,j and that of its projection 
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on the plane is cos(0j). Hence, to minimize the RI constant 
we should choose each cos(0j) as close as possible to 1, i.e. 

6*o = 9 1 = 02 = sin -1 J I that is actually the case reported 



in Figure [T] 

This choice clearly disregards the actual values of the 
probabilities pj and signal powers a 2 for j = 0, 1, 2 and may 
be suboptimal. 

To take these further information into account note that, 
since disturbances are introduced in acquiring the projections, 
they are independent from the plane. Hence, we may improve 
the SNR by choosing a plane that is able to rake a larger 
fraction of the signal power. We call this property rakeness 
and, in this case, to maximize it we have to maximize the 
power of the projection a 2 = po& 2 cos 2 (#o) +P1&1 cos 2 (6*i ) + 
P2&2 cos 2 (6» 2 ). 

With our assumption and by setting £j = cos 2 (6j) for j = 
0, 1,2 this amounts to maximizing a 2 = po<J 2 ^o + Pi<j\£,i + 
P2&2&> subject to the constraint on the 9j that becomes £ + 
£i+£ 2 = 2. Assuming thatpoc 2 . > Pi a i > P2&2> this criterion 
leads to £0 = £1 = 1 an d £2 = 0, i.e., a projections plane that 
coincides with the coordinate plane spanned by ao and a\. 

Clearly, the sheer maximization of the rakeness is not 
acceptable since any realization of a in which a 2 7^ would 
not be captured by the system or, in terms of the RI property, 
<5 = 1 since the a 2 axis belongs to the null-space of the 
projection operator. 

This toy case highlights that RI and rakeness may be 
suboptimal as a design criterion when considered alone and 
that improvements may be sought addressing the trade-off 
between RI enforcement and rakeness maximization. 

Such a trade off can be addressed both in a deterministic 
and in a statistical way. 

Pursuing the deterministic path, one may choose a pro- 
jecting plane like the one in Figure |2|-(a) that still allows 
signals along a 2 to have a non-zero projection but clearly 
favors directions with the largest expected power. Figure [2j- 
(b), that is analogous to Figure [T]-(b) for the new plane choice, 
shows that this is detrimental in terms of the RI constant 
since the length of the projection of the segment along a 2 
is substantially reduced. Yet, the lengths of the projections of 
the segments along ao and a\ are increased and since these 
are the occurrences carrying more power on the average, the 
overall average acquisition quality may be improved. 

The same improvement may be pursued in statistical terms 
by assuming that the projecting plane is chosen randomly at 
each measurement. In this case, the statistic of plane choices 
can be biased so that planes collecting larger energy are more 
probable, but planes allowing the acquisition of less important 
components are still possible. 

This second setting is particularly interesting since random 
projections are already employed to guarantee good RI prop- 
erties [6 | and the main aim of this contribution is to show that 
the trade-off between RI and rakeness can be addressed by 
proper design of the statistical distribution of the projecting 
directions. 

The rest of the paper is organized as follows. Section 
III will define the conversion architecture and lay down its 



mathematical model. Section [IV] introduces more formally the 
rakeness and its use as a design criterion. In doing this, to 
focus this exposition on application-oriented considerations, 
we accept that maximizing the energy of acquisitions is the 
right direction to go, thus postponing the statement of the 
formal chain of results starting from a mathematical definition 
of localization to a future contribution. This accepted, Section 
|V| describes a design path addressing the Rl/rakeness trade- 
off when the signals to acquire are localized in the frequency 
domain, that is by far the most common domain for signal 
analysis. Section VI expands that view to include a generic 



adaptive domain, that is able to reveal localization in a large 



class of signals. Section VII shows how the theory developed 
in the two previous Sections can be applied to the acquisition 
of signals like ECGs and small images of printed letters and 
numbers. Some conclusions are drawn at the end, while a 
couple of lengthy derivations are reported in the Appendix. 

III. System Definition 

We will concentrate on systems that perform AI conversion 
of sparse and localized signals by means of Random Modu- 
lation Pre-Integration (RMPI) (TJ. 

This scheme sketched in Figure [3] acts on signals of the kind 
a(t) where t is most usually time but may also be any other 
indexing variable. 

A "slice" of the signal a (say for —T/2 < t < T/2 for some 
T > 0) is processed by multiplying it by a waveform b(t) with 
a correspondingly sized support. The waveform b(t) is made 
by amplitude modulated pulses whose modulating symbols are 
chosen from a certain set. 

The most hardware friendly choices are rectangular pulses 
with binary ({0,1}) or antipodal ({— 1, +1}) symbols since 
multiplication can be implemented by a simple arrangement 
of switches that nicely embeds, for example, into switched- 
capacitor implementations Q. 

The resulting waveform is then integrated or low-pass 
filtered to obtain a single value that is converted into a 
digital word by conventional means. Note that this further step 
also nicely fits into a switched-capacitor implementation that 
naturally manages charge integration. 

Despite the fact that the rate (for time-indexed signals) or 
density (for generic signals) of the pulses may be very high 
and even larger than what a Nyquist-obeying acquisition would 
require, only the integrated values are actually converted. 

This multiply-and-integrate operation materializes the scalar 
product (a(t),b(t)) and can be performed M times (either 
serially or in parallel), each time considering a different 
waveform bj(t) {j = 0, . . . , M — 1) that is characterized by 
a set of modulating symbols drawn at random with a certain 
statistic. 

The resulting projections rrij — (a(t),bj(t)} for j = 
0, . . . , M — 1 can be aligned in a measurement vector m = 
(mo, . . . , rriM-i) T ■ Figure [3] exemplifies the signals and the 
operations entailed by the acquisition of rrij using antipodal 
PAM waveforms. 

For what concerns signal reconstruction, we assume that a 
is if-sparse, i.e., that there is a collection of N waveforms 
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Fig. 3. Block diagram of RMPI architecture: the signal to acquire is 
multiplied by the j-th antipodal PAM waveform and fed into an integrator 
whose output is sampled and quantized to produce the digital conversion of 
the j-th measurement 



To go further, we define the average rakeness p between 
any two processes a and j3 as 



,E 



a,0 



(3) 



where the constant n p is used to switch the meaning of p 
from "average energy of projections" (n p = 1) to "average 
power of projections" (e.g., n p = T -1 for signals observed in 

[-T/2,T/2]). 

It is worthwhile to highlight that p(a,/3) depends on how 
the second-order features of the two processes combine. In 
fact, we may expand the definition as 



Uj(t) for j = 0,. 
can be written as 



, N — 1 such that every realization of a 



a(t) 



N-l 

£ 

j=0 



CljUj(t) 



(1) 



for certain coefficients aj such that at most K < N of them 
can be non-zero at any time. 

Plugging ([T| into the definition of rtij we get rtij = 
Sfc^) 1 a k ( u k(t), bj(t)). By defining the vector a = 
(ao, . . . , aAr_i) T , the M x N projection matrix P = [P ■ k ] — 
[(uk(t),bj(t))], and the vector v_ = (vq, . . . , Vm-x) account- 
ing for the total noise affecting the projections, we have that 



Pa 



(2) 



is the reconstruction equality to be solved for the unknown 
a with the aid of its X-sparsity. In principle, this could be 
done by selecting, among all the vectors a satisfying Q, the 
one with the minimum number of non-zero entries. Since this 
is, in general, a problem subject to combinatorial explosion, 
many alternative theoretical and algorithmic methods have 
been developed allowing efficient and effective reconstructions 
[8 1 1 9 1 1 10]. Among all these possibilities, we will exploit the 
algorithm described in |9| in our experiments in Section VII. 

Note that y_ takes into account at least the intrinsic ther- 
mal noise affecting the analog processing of a(t) and the 
quantization noise due to digitalization. Since thermal noise 
is additive white and Gaussian (AWGN), its contribution to u 
is independent of the projecting waveforms bj(t) as long as 
they have constant energy. We assume that quantization noise 
is also approximately white and independent on the quantized 
input so that condition (ii) discussed in Section [TT] is satisfied. 

IV. Restricted Isometry and Rareness 

To cope with the noise in v, Rl-based design ifTTI tries to 
make the RI constant 5 of the projection operator as low as 
possible. This can be checked directly from the matrix P. In 
fact, since the projection is applied to isT-sparse vectors, we 
should consider each of the (^) matrices P' that are built 
selecting K of the N columns of P. If Ap, m and Ap, ax are 
respectively the minimum and maximum among the singular 
values 1121 of P' we have 



5 = max ■ 



{max[l-A£, in ,A|: ax -l]} 



p{a,0) = 



T/2 

-T/2 
T/2 

-T/2 



f a*(t)p(t)a{s)0*(s)dtds 

-T/2 J -T/2 
T/2 

E Q [a*(t)a(s)]Ep [f3(t)p* (s)]dtds 

-T/2 
T/2 

C a (t,s)C*p(t,s)dtds (4) 

-T/2 



where •* stands for complex conjugation and the two correla- 
tion functions C a and Cp have been implicitly defined. 

From the toy example in Section [n] we know that choosing 
the process b that maximizes the rakeness p(a, b) leads to 
good average SNR of the projections, but may destroy the 
RI property making the system insensitive to some signal 
components. 

To counter this over-tuning effect one may require that the 
process b is "random enough" to assign a non-zero probability 
to realizations that, despite being sub-optimal from the point of 
view of energy collection, allow the detection of components 
of the original signal that would be overlooked otherwise. 
Actually, this intuitive approach is fully supported by the 
existing results on the RI property. In fact, it is known J6l 
that if the matrix P is made of random independent entries, 
its RI constant is small with a substantially large probability. 

In general, enforcing the randomness of a process can be 
a subtle task since the very definition of what is random 
(entropic, algorithmically complex, etc.) can be extremely 
sophisticated and also dependent more on philosophical than 
technical consideration. 

Here, for simplicity's sake, we limit ourselves to en- 
ergy/power considerations and define a measure of the 
(non)randomness of a process as its self-rakeness, i.e., the 
average amount of energy/power of the projection of one 
of its realization onto another realization when the two are 
drawn independently. The rationale behind this quantification 
of randomness is that, if p(b, b) is high, then independent 
realizations of the process tend to align and thus to be 
substantially the same, implying a low "randomness" of the 
process itself H. 

This definition nicely fits into a mathematical formulation 
of the design path that increases the rakeness p(a, b) while 
leaving b random enough. In fact, given a certain sparse 
stochastic process a, we determine the stochastic process b 
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generating the projecting waveforms to employ in an RMPI 
architecture by solving the following optimization problem 



max;, p(a, b) 

(b,b)=e 
p(b, b) < re 



s.t 



(5) 



where e is the energy of the projection waveforms and r is a 
randomness-enforcing design parameter. 

Roughly speaking, solving |5]) will ensure that the result- 
ing waveforms will have constant energy (due to constraint 
< b,b >= e) paired with the ability of maximizing the average 
SNR of the projections (thanks to the capability of maxi- 
mizing the energy of the acquired samples since we impose 
that max;, p(a, b)) while maintaining the chance of detecting 
components of the original signal that carry smaller amounts 
of energy/power (thanks to the fact that each realization of the 
process b has p(b, b) < re, i.e., low autocorrelation and thus 
"large" randomness). 

In (|5j, the parameter e acts as a normalization factor, since 
if b' is the solution for e = e', then b" — ^Je" /e'b' is the 
solution of the same problem for e — e" . 

On the contrary, r is the parameter controlling the trade-off 
between the two design criteria we want to blend, i.e., RI and 
rakeness. Hence, different values of r lead to waveform with 
different final performance. 

Regrettably, though it is easy to accept that, thanks to their 
ability to maximize the energy of the samples, the resulting b 
may be able to increase the performance of the overall sensing 
system, the latter may rely (especially in the reconstruction 
part) on heavily non-linear and iterative operations that are 
difficult to model. For this reason, though feasible bounds 
for the parameter r can (and will) be derived theoretically 
in Section [V] and |VI[ the choice of its exact value is a matter 
of fine tuning of the global system, and it must be determined 
through numerical simulation. 



V. Localization in the frequency domain 

In this Section we specialize |5]l to the case in which 
the statistical features of a that cause the localization of its 
energy/power can be straightforwardly highlighted by Fourier 
analysis. 

We will concentrate on the time interval [— T/2, T/2] and 
set k p = T- 1 in Q. 

To express p in terms of the frequency-domain features of 
the processes a and j3 in it, let us assume that both of them 
are second-order stationary. 

Leveraging on this, we may define the single-argument 
correlation functions C a (s — t) = C a (t, s) and Cp(s — t) = 
Cp(t, s) whose Fourier transforms are nothing but the power 
spectra &(/) and /?(/) of the two processes. 

For p(a, 0) we obtain 



p(a,(3) 

I 



2T 



T-\P\ 



T+\P\ 



dqdp 



T poo poo 



C a (p)C^p)(l-^dp (6) 
«(/)/3*G?)e 27ri(/ ~ 9)p (l - |H dpdfdg 
a(fW(9) J* (l - dpdfdg 

&(f)0*(g)h T (f-g)dfdg (7) 



-T J — oo J — oo 
oo poo 



-oo «/ — oo 

OO pOO 



oo J —oo 



where 



sin 2 (7rT/) 



For simplicity's sake we may focus on the antipodal case 
in which the projection waveforms have a constant-modulus 
amplitude (±1), duration T, and thus automatically satisfy the 
constant energy constraint < b, b >= e in |5j with e = T, 
needed to make projection tuning possible. 

With this, the power spectrum of the projection waveforms 
can be designed by solving Q re-expressed in the frequency 
domain. To do so, use |7]i to rewrite p(a, b) and p(b, b) in |5]) 
and consider 



oo poo 



max 

b j — oo «/ — oo 



oo />oo 



&(f)b(g)h T (f - g)dfdg 
b(f)b(g)h T (f-g)dfdg<rT 



oo J —oo 



Kf) > o 



(8) 



s.t. 



KfW = i 



Kf) = K-f) 

where the last three constraints encode the fact that b must be 
a power spectrum of a unit-power, real signal. 

Once that r is fixed, ((HJ can be solved by assuming that 
a concentrates its power in the frequency interval \—B, B] 
and applying some kind of finite-elements methods, i.e., ap- 
proximating all the entailed functions with linear combinations 
of basic function elements on which the integrals can be 
computed at least numerically. 

As an example, select a frequency interval [— B, B] and par- 
tition it 2n + l subintervals of equal length A/ = 2B/(2n+l) 
Fj = [j - A//2,j + A//2] for j = -n,...,n. Assume 
now that &(/) is constant in each Fj and define Xj(f) as 
the indicator fucntion of Fj, i.e. Xjif) = 1 if / S -Fj and 
otherwise. We have b(f) = X}"__ ra bjXj (/) f° r certain 
coefficients 

This can be substituted into |8]l to obtain 
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max 

b_„,...,6„ 



s.t. 



with 



E w i b i 

j=—n 

n n 

E E WkW jlk < tt 

j——n k——n 

bj > j — — n, . . . , n 

n 

Af e h = i 

3=—n 

by = b-j j = -n, . . . , n 



a(f)h T (f - g)dfdg 



h T (f-g)dfdg 

F k 



(9) 



This leaves us with the vector of 2n+l unknown coefficients 
&_„, . . . , b n that must determined by solving an optimization 
problem characterized by a linear objective function and few 
linear and quadratic constraints. Plenty of numerical methods 
exist for solving such problems even for large number of basic- 
elements and thus for extremely effective approximations 
(commercial products such as MATLAB or CPLEX provide 
full support for large-scale version of these problems). 

Once that the optimum b(f) has been computed, one may 
resort to known methods to generate an antipodal process with 
such a spectrum exploiting a linear probability feedback (LPF) 
|f]~3l lfl4l |fl5l . Slices of length T of this process can be used 
as projection waveforms in an RMPI architecture for the CS 
of the original a. 

Note that, even if this is needed to arrive at a final working 
system, the core of rakness-based design concerns the solution 
of |5) for frequency-localized signals to obtain the best spectral 
profile of the projecting waveforms, independently of their 
physical realization. How such a spectral profile can be ob- 
tained using antipodal PAMs is an implementation-dependent 
choice, which allows to realize an hardware system for sparse 
and localized signal acquisition which does not require analog 
multipliers 0. 

As far as the range in which r should vary to administer the 
trade off between RI and rakeness, note that, since p(b, b) is a 
measure of (non)randomness, it must be minimum when the 
process b is white in its bandwidth, i.e., when b(f) — 1/(2B) 
for / € [— B, B] and otherwise. 

Plugging this into (|7]) and defining c = BT one gets 



ing function of c with lim c _ i . 00 r mm c = 1/2. Hence, we may 
safely use such an asymptotic value to set l/(2c) as a suitable 
lower bound for r in any practical conditions. 

Again, from the meaning of p(b, b) we got that it is 
maximum when the waveforms produced by the process are 
constant. This implies Cfc(r) = 1 that can be plugged into |6]) 
to obtain 



r < r max = 



1 



T 



dp = 1 



Overall, the tuning of the overall system will optimize 
performance by choosing r E \^~>M- 



VI. Localization in a generic domain 

Slices of second-order stationary processes (that enjoy a 
simple and well-studied characterization in the frequency 
domain) do not exhaust the set of signals that we may want 
to acquire. 

To cope with more general cases assume to work in 
normalized conditions such that both the waveforms to be 
acquired and the projection waveforms have unit energy, i.e., 
f\ \a(t)\ 2 dt = JJt \b(t)\ 2 dt = 1, where the latter constrain 
sets e = 1 in |5]). 

When we comply with this assumption (possibly by scaling 
the original signals), if C x represents either C a or we have 
that 

• C x is Hermitian, i.e., C x (t,s) = C*(s,t); 

• C x is positive semidefinite, i.e., for any 
integrable function the quadratic 

form J% f% C (t)C x (t, s)Z(s)dtds 

2 2 

21 



E 



l\ x(i)£(«)df 



yields a non-negative result; 



has a unit trace, i.e. 



C x (t,t)dt 



E[|cc(£)| 2 ]di = E 



\x(t)fdt 



From this, we know (see e.g. Ifl6l ) that two sequences of 
orthonormal functions 6o(t), 9i(t), . . . and (f>o(t), (/)%(£),... 
exist, along with the sequences of real non-negative numbers 
Mo > A*i > • • • an d Ao > Ai > . . . such that ^°^ Mi = 

J2T=0 X 3 = 1 and 



r > r min = 

Ci(47rc) + 47rcSi(4-7rc) - log(47rc) + cos(47rc) - 7 - 1 

47T 2 C 2 

where 7 is the Euler's constant and Ci and Si are respectively 
the cos-integral and sin-integral functions. 

The quantity r mm c is a monotonically and rapidly increas- 



c a (t, s ) = 5>^;(i)^-( s ) 



3=0 



c b (t,s) = E^-W'M* 
3=0 



(10) 

(11) 



By substituting the generalized spectral expansions for the 
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two correlation functions ( |T0] > and ( fTTj ) into Q one gets 



oo oo 



P( a > h ) = E E X i^j,k 

j=0 k=0 
oo 

p(6,6) = 
where the real and nonnegative numbers 



^j,k 



are the squared modulus of the projections of each <frj on every 
Ok (and viceversa). 

The orthonormality of the Ok guarantees that the sum of 
the squared modulus of the projections of <f>j must equal the 
squared length of <fij itself and thus, since <f>j is normal, that 
jy^Lo ^i. fe = Conversely, from the fact that the ijjj are 
orthonormal we have also X^fcLo = 1- 

Hence, the optimization problem |5]l can be rewritten in 
totally generic terms as 



oo 



max max / y / y ^jMfc^j.fc 



j=0 fc=0 



oo 



S.t. 



3=0 

OO 

E^ 2 

3=0 



^3,k 
oo 

E E j-k 

3=0 

OO 

E E i> k 

k=0 



> Vj 

= 1 

< r 

> Vj,fc 
= 1 Vfe 

= 1 Vj 



(12) 



Note that the two max operators address separately the 
problem of finding an optimal basis (max^) and then the 
optimal energy distribution over that basis (max^). 

As far as the range of r is concerned, assume to know that 
J is an integer such that Xj = for j > J. It can be easily 
seen that max^^Tg 1 X 2 j subject to the constraints Xj > and 
^2j=q Xj — 1 is 1 and is attained when Ao = 1 and Xj = for 
j > 0. It is also easy to see that min J2j=o A| subject to the 
constraints Xj > and J2j=o A^ = 1 is 1/ J and is attained 
when X 3 = 1/J for j = 0, . . . , J - 1. Hence, r € [1/J, 1]. 

In particular, the lower bound r > 1/J rewritten as rJ > 1 
can be read as a general rule of thumb, i.e., the more random 
the process that generates the projection waveforms, the larger 
the number of non-zero eigenvalues in the spectral expansion 
of its correlation function. 

The solution of ( fl2] i is derived in the Appendix and depends 



on the two partial sums 

Ei(J) 

s 2 (J) 

to obtain 



j-i 

3=0 
.7-1 

3=0 



Xj = Xj(J) = 



J/ij - Si(J) 
's 2 (J)-iEf(J) 



(13) 
(14) 

(15) 
(16) 



which hold for j = 0, 1, . . . , J — 1 where J is defined by 

J = max{j | Xj-iU) >0} (17) 
By definition, all the eigenvalues Xj for j > J are null. 

A. Finite dimensional signals 

The special case in which the signal to be acquired can 
be written as a linear combination of known waveforms by 
means of random coefficients is, for us, extremely interesting 
and deserves some further discussion. 

Let us assume that ([TJ holds for orthonormal Uj (j = 
0, . . . , N — 1) and let us compute 



JV-l N-l 

a *j a k]Uj(t)u k (s) 

j=0 fc=0 



(18) 



The correlation matrix A 



[A j>k ] = [E[a*a k }] is Her- 
mitian and positive semidefinite, hence it can be written as 
A = QMQ? where ^ stands for transposition and conjugation, 
M is a diagonal matrix with real non-negative diagonal entries, 
and Q is an orthonormal matrix whose columns are the 
eigenvectors of A. 

With this, we may rewrite (fl8]l as 



N-l N-l JV-l 



C a (t,s) 



u*j(s)u k (t) 



j=0 k=0 1=0 
JV-l JV-l JV-l 

;=0 .3=0 k-0 

Hence, we may express C a (t,s) in the form needed for 
writing flO} and thus the solution of ( fl2] > by simply setting 

°3 = Efc=o Q*k,j u k and Mj = Mjj for j = 0, . . . , N - 1. 

This straightforward derivation clarifies that, when we have 
identified sparseness along a certain signal basis, the statistic 
of the coefficients gives us hints on the basis that may be used 
to highlight localization. Along this other basis, localization 
itself is nothing but the difference between the lower-index, 
largest eigenvalues /io> Mi> • • • an d me others. 

A bridge is also built between the general treatment of 
rakeness in this Section and the frequency-domain analysis 
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of the Section [V] In fact, if a is substantially bandlimited in 
the frequency interval [—B,B] and is considered in the time 
interval [— T/2, T/2] its realizations may be well expressed 
as a linear combination of waveforms that are the truncated 
version of prolate spheroidal wave functions I:il7l fl8l . It is 
known that, if c = BT then N — 2c functions are enough to 
achieve an approximation quality that dramatically increases as 
c — > oo. Hence, the solution of ( p~2] > will feature J = N = 2c 
for values of r G [1/J, 1] = [y-, 1]. 

From an operative point of view, whatever analysis allows 
us to obtain the generalized spectral expansion of C a as in 
( fT0l >, we may use ( fl3] l, ( p~6] >, ( fTTj i and ( fTT) to compute the 
correlation function Cb of the process generating the projection 
waveforms. 

To fit this Cb into an actual RMPI architecture, we must 
generate a binary or antipodal PAM signal with such a non- 
stationary correlation. The details of the mechanism allowing 
this are far beyond the scope of this paper and will be the 
topic of a future communication. 

It is here enough to say that, if the number of symbols 
S in each waveform is limited to few tens (say S < 100), 
we are able, depending on C&, to automatically determine 
two sets of cardinality s = S(S + l)/2: the first set 
{zq, zi, . . . , z s -i} contains sequences of modulating symbols, 
while the second set {£o, Ci, ■ ■ ■ , Cs-i} contains probabilities, 
so that = 1- 

These two sets are such that, if each time a projection 
waveform is needed, the modulating symbols in Zj are used 
with probability Q, then the resulting process has the desired 
correlation. 

In any case, let us stress that, as noted before for frequency- 
localized signals, the core of rakeness-based design concerns 
the solution of (|5), which is here described for generically lo- 
calized signals. Once that the correlation of the best projection 
waveforms is determined, their actual realization depends on 
implementation assumption that may vary from application to 
application. 

VII. Sample applications 

In this Section we introduce rakeness as a design criterion 
to optimize the performance of two acquisition systems, one 
that deals with Electro Cardio Graphic (ECG) signals, which 
can be easily modeled in the frequency domain as we did in 
Section [V] and the other that must be described relying on the 
generalized spectral expansions in Section |VI] since its target 
signals are images. 

Despite the fact that the two scenarios are different, the 
path we follow in designing an acquisition system based on 
CS system is the same and can be summarized in few steps: 

i) identify the basis with respect to which the signal to 
acquire is sparse; 

ii) identify the basis with respect to which the signal is 
localized; 

iii) solve |5]l for a number of possible values r in its range; 

iv) for each value of r, implement an RMPI architecture 
exploiting the sparsity revealed in and in which the 
projection waveforms are as close as possible to the 
optimal ones; 



v) perform Monte-Carlo simulations to evaluate the result- 
ing systems and select the best performing one. 

Note that, in the classical design flow of a CS system, |TJ> is 
a prerequisite while [Iv} is tackled once assuming that the pro- 
jection waveform are random PAM signals with independently 
and identically distributed ("i.i.d." from now on) symbols. This 
is what will be taken as the reference case to quantitatively 
assess the improvements due to rakeness-based design. 

In all cases, the performance index is the average recon- 
struction SNR (ARSNR), i.e. the average ratio between the 
energy of the original signal over the energy of the difference 
between the original signal and the reconstructed one. ARSNR 
values are always plotted at the center of an interval accounting 
for the variances of the corresponding reconstruction SNRs. 

Note also that the implementation constraints (e.g., the 
restriction of projection waveform to PAM profiles with an- 



are nothing but an 



in 


iv 


in 


V 



Finally, one may observe that steps 
elementary line-search for the best possible value of r. As a 
matter of fact, the values of r for which a definite improvement 
can be obtained are easily identified by means of a very small 
numbers of trials. 



A. Acquisition of ECGs 

The ECG time shape represents the voltage between two dif- 
ferent electrodes placed on the body at two specific positions. 
It records the electrical field produced by the myocardium, 
i.e., for each heart beat it cyclically reports the successive 
atrial depolarization/repolarization and ventricular depolariza- 
tion/repolarization. 

The application of CS techniques to ECG acquisition has 
been the topic of recent contributions 1 19ll20ll 2T1ll22l aimed 
to either reducing the amount of data needed to represent the 
signal in mobile applications or to achieve an high compres- 
sion ratio in data storage systems. 

In order to demonstrate the effectiveness of rakeness-based 
design for CS of ECGs, we need a broad collection of 
realistic realizations. To achieve this goal, we used a synthetic 
generator of ECGs, thoroughly discussed in ll23ll that provides 
signals not corrupted by noise to which we add white Gaussian 
noise with suitable power. The amount of noise is chosen 
so that the considered environment is realistic, but it can be 
arbitrary from the point of view of rakeness-based design that 
is independent of the noise level. 

The generator core is expressed by the following set of three 
coupled ordinary differentially equations ll23ll 



Xi 
±2 

%3 



U\Xi — U>2X 2 
UJ1X1 + UJ2X2 



i£{P,Q,R,S,T} 



7,6j exp 



2v? 



(x 3 - x 3 ) 



(19) 



Each heart beat is represented by a complete revolution on an 
attracting limit cycle in the (xi, £2) plane. The shape of ECG 
signal is obtained introducing five attractors/repellors points 
in the x 3 direction in correspondence to the peaks and valleys 
that characterize the time shape of the signal and which are 
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TABLE I 

Parameters bounds used in the ECG generator. 



Intrinsic SNR = 17 dB 



Index 


P 


Q 


R 


S 


T 


0* 


-75;-65 


-20;-5 


-5;5 


10;20 


95;105 


7» 


1;1.4 


-5.2;-4.8 


27;33 


-7.7;-7.3 


0.5;1 




0.05;0.45 


-0.1;0.3 


-0.1;0.3 


-0.1;0.3 


0.2;0.6 



<D 12 

Q 
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g 
o 



H 




ECG signals 


■ ■! 

• i 




- - - high r 


• • 




optimum r 


■ ■ 




™ ™ ■ low r 


• ■ 






• ■ 






■ ■ 






■ ■ 






• ■ 






■ ■ 






■ ■ 


























s; .yg* 





0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

Normalized Frequency 



Fig. 4. Average spectra of real ECG signals (gray area) and of the sampling 
PAM sequences corresponding to the optimum (solid line) as well as an high 
(dashed line) and a low value (dash-dotted line) of r. 



conventionally labeled by P, Q, R, S and T; furthermore a; 3 in 
( fT9] > represents the mean value of the generated ECG. 

In order to mimic the behavior of ECGs in patients affected 
by the most studied cardiac illness, the parameters ji, Vi 
and Qi, i G {P,Q,R,S,T}, characterizing each considered 
signal are taken from a set of random variables uniformly 
distributed within the bounds reported in Table [I] In addition, 
we randomly set the heart rate between 50Hz and 100Hz by 
property adjusting loi and oj 2 - 

Though CS methods are classically developed for sparse 
representation with respect to signal bases, they have a 
straightforward generalizzation to sparse representation with 
respect to dictionaries, i.e., redundant collections of non- 
indipendent waveforms I24l ll25ll . 

This is, in fact, the case of ECGs, for which a dictionary 
made of Gabor atoms 



J 

Qs.u.v .w\J') 

V s 



-0 



cos(vt + w) 



can be used (26) (27). 

In our experiment, a total of 507 atoms are used corre- 
sponding to different quadruple of parameters (s, u, v, w). This 
collection of Gabor atoms is obtained using a greedy algorithm 
able to extract a limited number of functions from a broader set 
[27 1. With respect to this dictionary the sparse representation 
of a typical ECG heartbeat waveform requires about 14 non- 
zero coefficients. 

Furthermore, in our simulations, T = Is within which the 
signal is sampled N = 256 times (a common choice for ECG 
equipments). 

To apply the results in Section [V] we first compute the aver- 
age of the power spectral densities of 1000 ECGs generated as 
presented above to obtain a(f), the input of the optimization 



CO 

g 13 





O ave ± std 
- - -i.i.d. sequences 

^^^h lnrpli7Prl QpniipnppQ 

^^^^^ lUwl I^CU oCljUCI IVCO 


5 -s-.. 










•••••i - 



7 

N/M 



Fig. 5. Average value of the reconstructed SNR (ARSNR) as a function 
of the signal compression ratio N/M between the number of Nyquist samples 
and of CS measures. The dashed line refers to i.i.d. sampling waveforms and 
the solid line to rakness-optimized ones. 



problem ([HJ. The shape of a(f) is the gray profile shown in 
Figure [4] Next, we find the optimum r as described at the 
beginning of Section VII Figure [4] shows the optimum profile 
for the best value r = 0.038 (solid line) as well as for a smaller 
value (dash-dotted line) and for a larger value (dashed line). 

Finally the LPF generator mentioned in Section[V]is used to 
produce the antipodal sequences with the optimized spectral 
profiles. These sequences are used to take M measurements in 
a time window of length T, to which we add white Gaussian 
noise to construct the measurement vector m according to Q. 

To determine the performance of rakeness based design, we 
consider a test set of 2000 ECG signals different from those 
employed for determining the average spectrum of Figure 
[4] These signals are acquired by projecting them both on 
localized antipodal sequences and on i.i.d. antipodal sequences 
(classically employed in CS-based methods and our reference 
case). The resulting ARSNRs are shown in Figure [5] as a 
function of the ratio between the intrinsic dimension of the 
signal and the number of CS measures. In both cases the 
intrinsic SNR is equal to 17dB. 

As it can be noticed, rakeness-based design allows to 
achieve an improvement of at least 3.5dB in ARSNR with 
respect to the i.i.d. case, and even yields denoising (i.e. and 
ARSNR larger than the intrinsic SNR) for small compression 
ratio values. To give a visual representation of the improve- 
ment, Figure [6] reports, for M = 32, a comparison between 
an ECG signal and the reconstructed one for the i.i.d. (a) and 
rakeness-based (b) case. Direct visual inspection is enough to 
confirm the superiority of our approach. 



B. Acquisition of small images 

The signal to acquire is a 24 x 24-pixel image, each pixel 
value ranging from (black) to 1 (white), which represents a 
small white printed number or letter on a black background 
with a gray-level dithering to make the curves smoother to 
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Fig. 6. Original (solid line) and reconstructed (dashed line) ECG when 
i.i.d sampling waveforms are used (plot (a)) and when rakeness-optimized 
sequence are exploited (plot (b)). In both cases N=256 and M=32 and the 
intrinsic SNR=17dB. 



the human eye. Number and letters are randomly rotated and 
offset from the center of the image but never clipped. 

Although due to random rotations and offsets almost all 
pixels have a non-vanishing probability of being non-zero, a 
typical image contains only about 85 bright pixels, so that 
can be considered sparse in the base of 2-dim discrete delta 
functions that evaluate to 1 at a single pixel position and zero 
elsewhere. We may thus think of acquiring them using a RMPI 
architecture that projects along 24 x 24 antipodal random grids 
to obtain measurements that are enough to reconstruct the 
image but whose number M is much less than the number 
N = 24 x 24 = 576 of the original pixels. 

To simplify the design phase, the generation of the random 
grids is done by adjoining 4x4 subgrids each with 6x6 
antipodal values whose statistic is optimized by solving (12) . 

To allow calculations, the values in each subgrid are rear- 
ranged into a 36-dimensional vector as schematically reported 
in Figure [7] that also highlights the subgrid on which we will 
focus in the following. 

In that region, and due to the vector rearrangement, we may 
list the modulating symbols of the projection waveform b with 
bj for j = 0, . . . , 35. The same can be done for the incoming 
signal a when it is expressed along the basis of 2-dim discrete 



■ 

■ 
■ 


■■■ 


■ 


■i 





Fig. 7. A sample image, its partition and rearrangement into a vector 
containing the value of each pixel. 



delta's with coefficients that we may indicate with a,j for j = 
0, ...,35. 

If a = ( ap, ■ ■ ■ , a35) T , we may follow the development of 

EfaoH 



subsection 



VI-A 



to estimate the 36 x 36 matrix A 

by empirical averaging over a training set of 2080 randomly 
generated images. 

The resulting matrix is reported in graphic form in Figure 
jSJ-Ca) where, for each pair of indexes j, k = 0, . . . , 35, a point 
is laid down whose brightness is proportional to the values of 



E[a,-a fe l. 



The eigenvalues (1q , . . . , /i 35 of that A are reported as the 
light bars in Figure [8]-(c). By exploiting ( [16) and ( [T7| ) for 
r = 0.047 we get J = 36 and the eigenvalues Ao, ■ ■ ■ , A35 
reported as dark bars in figure [8]-(c). 

From the eigenvectors of A and these new eigenvalues, 
we may construct the correlation matrix of the values in this 
projection subgrid. Since the subgrid contains 6 x 6 = 36 
values, its correlation matrix 13 has dimensions 36 x 36 and 
is reported in Figure [8]-(b) in a graphical form adopting the 
same convention used to represent the values of A. 

Once that B is known, we are able to select a collection of 
36 x 35/2 = 630 antipodal subgrids Zq, ... , Z629 with attached 
probabilities (q, • • • j C629 so that, if each time a projection is 
needed we use the subgrid Zj with probability Q, the overall 
process features exactly that correlation matrix. 

The same design process is repeated for each of the 4 central 
6x6 regions in the image while the 12 outer subgrids are built 
from independent and uniformly distributed antipodal sym- 
bols. All the subgrids are finally compounded in a complete 
24 x 24 projection grid. 

As a comparison case, projections are also taken by using 
i.i.d. symbols for all the elements of the projection grid. 

In both cases, noise is added to the projections before they 
take their place in the vector m of M measurements according 
to |2]), and reconstruction is performed using the algorithm 
reported in [9]. 

Figure [9] reports the ARSNR (over 3000 trials) of the 
reconstructed images and compares the performance of an 
RMPI based on rakeness-optimized projections and on i.i.d. 
projections for different values of the compression ratio N/M. 
In both cases the intrinsic SNR is 17dB. 

It is evident from Figure [9] that, even if it is exploited only 
in the central portion of the images, rakeness-based design 
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Fig. 8. Correlation matrix of the pixels in one of the four central regions 
(a) Correlation matrix of the optimal projection processs (b) Eigenvalues of 
the above correlation matrices (c). 
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Fig. 9. Quality of the reconstructed images when rakeness-optimized or i.i.d. 
projection grids are used in an RMPI architecture for different compression 
ratios. 



leads to non-negligible improvement of at least ldB. 

A qualitative appreciation of such an improvement can be 



obtained from Figure 10 in which 5 images (a) are acquired 
and reconstructed by means of M = 115 rakeness-optimized 
projections (b) or by the same number of i.i.d. projections 
(c). Reconstruction artifacts are visibly reduced by adopting 
rakeness-based design. 




Fig. 10. Sample images (a) and their reconstruction based on rakeness- 
optimized projection grids (b) or on i.i.d. projection grids (c). 



VIII. Conclusion 

Compressive sensing exploits the fact that, when looked at 
in the right domain, the information content of a signal can 
be much less than what appears when we look at it in time or 
frequency (i.e., the signal is sparse). 

Acquisition schemes that exploit sparsity may lead to con- 
siderable advantages in terms of sensing system design since, 
for example, if the information content is much less than the 
signal bandwidth, sub-Nyquist sampling can be employed. 

To all this we add the consideration that, in a possibly dif- 
ferent domain, the energy of the signals may be not uniformly 
distributed (i.e., the signal is localized) and when noise is 
present, it is convenient to adapt the system to "rake" as much 
signal energy as possible. 

By itself, this is not a novel concept since it appears, 
for example, in matched filters and rake receivers used in 
telecommunications. Yet, in our context, the efforts to collect 
the energy of the signal must be balanced with the guarantee 
that all details of its underlying structure can be captured 
when immersed in noise. This brings us to a trade-off that 
we propose to address in statistical terms by means of an 
optimization problem: maximize the "rakeness" while obeying 
to a constraint ensuring that the measurements are random 
enough to capture all signal details. 

The paper develops the formal definition of such problem as 
well as its solution for stationary signals whose localization 
can be highlighted in the frequency domain, and for more 
generic non-stationary signals whose localization is more 
evident in suitably defined domains. 

The applicability of both techniques is demonstrated by 
sample applications to the acquisition of ECG tracks and small 
letter images. 

IX. Appendix 

A. Solution of ( |12| > 

The first subproblem max^ can be solved leveraging on the 
fact that it is a linear problems with linear constraints. Since, 
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in principle, it may involve an infinite number of variables we 
should proceed by steps. 

Let P n be the optimization problem 



we have 



max 



17-1 n- 1 

\ " \ " 



j=a fc=o 

Sj t k > o yj,k 

s.t. E^lo^i.fe = 1 Vfc 

so that Pqo is the basis finding subproblem in ( |12) , 

Since all the series involved in the definition of Poo are 
convergent, we have that, independently of 3^ k, 



7*— 1 77— 1 

lim EE 

j=0 fc=0 



OO OO 



Xj^k^j.k 2^ 2^ A,//;,- ,./, 

j=0 fc=0 



Moreover, since all the summands are positive, the limit is 
from below. 

Let us now assume to have solved Poo yielding a value 
a (Poo) corresponding to a certain optimal choice 3°° fe . 

Given any e > there is a n such that for any n> n 

71—1 71 — 1 

< a(Poo) -EE V*S~fc < e 
Yet, by solving P n we get a solution u(P n ) such that 



71— 1 71— 1 

E E n 3t~*-^i 

j=0 fc=0 



X jf x k 3^ k < a(P n ) < ^(Poo) 



where the last inequality holds since every feasible configura- 
tion for P n is also a feasible configuration for Poo- 

Altogether we get that for any n> n 



that is 



< a(Poo) - ct(P„) < e 
lim cr(P„) - (j(Poo) 



from below. 

From this we know that, if the solutions S" fc of P n have a 
limit, such a limit yields ct(Poo)- 

To study the solutions of P„ we may first recall that the 
polytope 



Sj.fc > j, fc = 0, . . . , n - 1 

YTjZo 3j,fe = 1 = 0, . . . , n - 1 
J2lZo E J,k = 1 i = 0, . . . , n - 1 



is the one characterizing the so-called "assignment" problems 
[28 1 and is well known [29] to have vertices for 3^ for j, k — 
0, . . . , n — 1 equal to a permutation matrix. Hence, let £ : 
{0, 1, . . . , n — 1} M> {0, 1, . . . , n — 1} be the bijection such 
that 

Jl iffc = £(j) 
otherwise 



71-1 



3=0 

for some optimally chosen £. 

Actually, we may prove that such an optimal £ is the 
identity. We do it by induction. 

For n = 2 there are only two permutations corresponding 
to the two candidate solutions a' = Ao/io + ^lfa and rr" = 
Ao^ti + Ai/io- Yet, from the sorting of the Xj and of the /ij 
we have a' — a" = (A — Ai)(/xq — Mi) ^ 0- 

This confirms that the optimum solution is the one corre- 
sponding to £(j) = j for j = 0, 1. 

Assume now that this is true for n up to a certain n and 
that we have solved P ft +i by means of a permutation £. 

If £(0) = J> then er(P,- i+ i) = AoMj+cr'. Yet, cr' must be 
the value of the solution of a problem with n terms Ai, . . . , A fi 
and fix, ... , fij-i, Mj+1j ■■-.ft- Since we assumed to know 
how problems with n terms are solved we know that 



3 

E 



J=J+1 



It is now easy to see that the value Xofij + a' of the alleged 
solution is actually smaller than YZj=o fa- 
in fact 



E ^3 ^3 - ^Ofa - E ^3'fa-l 



3=0 



3=1 



E X 3fa 
3=3+1 



3 



Ao(w> - Mj) - 2^ ^(/ij-i - Mi) 

3 

fa) -^2 X 3(fa -fa-l) 



3 

A o E(w- 



3 = 1 



3 = 1 



3 

E 

.7=1 



(Ao-Aj)(Mi-i -Mj) >0 



Hence, the optimal permutation must feature £(0) = 0. This 
reduces the solution of P fi +i to the solution of P ft that we 
already know to be = j for j = 1, . . . , n — 1. 

In the light of this, every P„ has a solution corresponding 
t0 £(i) = J f° r J — 0, • • ■ , n — 1 and the solution of P^ is 

3j,fc = 

GO 

er(Poo) = E^'W 

J=0 

This solves the basis-selection problem and yields ( fT5| ). 
The original ( fT2] i now becomes 



max a Aj/ij 
j=o 
A 



> Vj 
^•t- E, x ., A., = 1 

< r 



.7=0 A J 



(20) 



'3=U ' J 

Since the Xj are non-negative and sorted in non-increasing 
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order we have that the set of indexes such that Aj > must be 
of the kind {0, 1, . . . , J — 1} for some integer J > 0. We also 
know that, to allow J2j = Q Aj = 1 and Y^,j=o Aj = r to hold 
simultaneously we must have r > 1/J and thus J > 1/r. 
Hence, for a given J > 1/r our problem can be recast into 



J-l 

max> 



Such equations can be solved for £' and £" and the resulting 
values substituted in ( |22"1 ) to yield ( fTo) . 



s.t. 



/ , - -3 Mj 

3=0 

Ej=o A-. 



> 3 
= 1 



0,...,J-1 



E^A| < 



(21) 



B. Real and nonnegative values in ( |16| l 

The denominator within the square root is positive whenever 
r>l/J. 

To show that the corresponding numerator is also non- 
negative write 

J£ 2 (J) - £ 2 (J) = 

J-l ,7-1.7-1 
3=0 3=0 k=0 



J-l 



J 



Note that the feasibility set of ( f2"T] i for a certain J = J 
contains points that are arbitrarily close to those of the 
feasibility set of pi) for any J < J. Hence, to maximize 
the rakeness we should try to have J as large as possible. 

To determine the J leading to maximum rakeness note first 
that, if we drop the randomness constraint — r > tne 

relaxed problem has the trivial solution A = 1 and Aj = for decreasing and such that J2jZo (j = 0. Hence, there is a f 



J-i 

3=0 



k=0 



Let now Q = /Zj — j J2k=o Mfc- We have that the Q are 



j > 0. Such a solution is not feasible for the original problem 
since Ylj=o A| = 1 > r, hence the corresponding optimum 
must be attained when the randomness constraint is active, i.e. 
for EjZo A ■ = r. 

The Karush-Kuhn-Tucker conditions for ( pO) with the in- 
equality constraint substituted by the equality constraint are 



+ e + i"x. 



S)j=0 ^3 

Eoo ,2 
3=0 A j 

l 3 

l 3 A l 



Vj 
> Vj 

= 1 



> Vj 
= Vj 



where £' is the Lagrange multiplier corresponding to 
Y^jLo Aj = 1, I" is the Lagrange multiplier corresponding 
to £°lo\ 2 = r, and ^ 
corresponding to Aj > which must hold Vj 



; , and l"J are the Lagrange multipliers 



Since for Aj > the constraint I j A j 



sets 



know that 



A 3 - ^r- 



we 



(22) 



for j = 0, . . . , J — 1. 

Since the sequences Aj and fij are both decreasing, we must 
have I" < and thus I' > -fij for j = 0, 1, . . . , J - 1, i.e., 
£' > — fJ,j-i. Hence, 

J = max{j | I > — fjtj-i} 

To see how the two parameters £' and £" depend on J note 
that they should satisfy the simultaneous equations 



Z^j=0 -£" 
2-^3=0 \ I" J 

i.e., exploiting ( [T3| ) and (14) , 

Si(J) + Jf 



1 

r 

= -f 



S 2 (J)+2Si(J)f + Jf 2 = /i" 2 r 



such that EjLo 1 C3 = E/=/'(-&) > °- Hence ' 

JE 2 (J)-S 2 (J) = 

i'-i .7-1 

_3=0 3=3' 

j'-l J-l 

J Mi'-i£C?-A*j'XJ(-Cj) 
3=0 3=3' 

i'-i 

J(Mj'-i -Mj') X C > 

3=0 
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